## pval_cutoff: 0.05
## lfc_cutoff: 1
## low_counts_cutoff: 10
General statistics
# Number of samples
length(counts_data)
## [1] 6
# Number of genes
nrow(counts_data)
## [1] 55487
# Total counts
colSums(counts_data)
## SRR13535276 SRR13535278 SRR13535280 SRR13535300 SRR13535302 SRR13535304
## 3107284 2321609 3701956 2491487 1580539 1861995

Create DDS objects
# Create DESeqDataSet object
dds <- get_DESeqDataSet_obj(counts_data, ~ treatment)
## [1] TRUE
## [1] TRUE
## [1] "DESeqDataSet object of length 55487 with 0 metadata columns"
## [1] "DESeqDataSet object of length 14648 with 0 metadata columns"
colData(dds)
## DataFrame with 6 rows and 25 columns
## Assay Type AvgSpotLen Bases BioProject BioSample Bytes Center Name Consent DATASTORE filetype DATASTORE provider DATASTORE region Experiment treatment GEO_Accession (exp) Instrument LibraryLayout LibrarySelection LibrarySource Organism Platform label ReleaseDate Sample Name source_name SRA Study
## <character> <numeric> <numeric> <character> <character> <numeric> <character> <character> <character> <character> <character> <character> <factor> <character> <character> <character> <character> <character> <character> <character> <factor> <POSIXct> <character> <character> <character>
## SRR13535276 RNA-Seq 300 8225466000 PRJNA694971 SAMN17588686 3252113587 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943360 A GSM5043430 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA in space without gravity 2021-09-09 GSM5043430 C2C12 proliferating .. SRP303354
## SRR13535278 RNA-Seq 300 9203426700 PRJNA694971 SAMN17588684 3619152333 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943362 A GSM5043433 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA in space without gravity 2021-09-09 GSM5043433 C2C12 proliferating .. SRP303354
## SRR13535280 RNA-Seq 300 9323939700 PRJNA694971 SAMN17588682 3735905901 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943364 A GSM5043436 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA in space without gravity 2021-09-09 GSM5043436 C2C12 proliferating .. SRP303354
## SRR13535300 RNA-Seq 300 12820015200 PRJNA694971 SAMN17587361 5047533646 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943384 E GSM5043471 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA on land 2021-09-09 GSM5043471 C2C12 proliferating .. SRP303354
## SRR13535302 RNA-Seq 300 12499917600 PRJNA694971 SAMN17587359 4941074444 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943386 E GSM5043475 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA on land 2021-09-09 GSM5043475 C2C12 proliferating .. SRP303354
## SRR13535304 RNA-Seq 300 7150086300 PRJNA694971 SAMN17587357 2845819297 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943388 E GSM5043478 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA on land 2021-09-09 GSM5043478 C2C12 proliferating .. SRP303354
Sample-to-sample comparisons
# Transform data (blinded rlog)
rld <- get_transformed_data(dds)
PCA plot
pca <- rld$pca
pca_df <- cbind(as.data.frame(colData(dds)) %>% rownames_to_column(var = 'name'), pca$x)
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 40.9653 35.5987 26.4305 19.10170 17.10922 6.531e-14
## Proportion of Variance 0.3901 0.2946 0.1624 0.08482 0.06805 0.000e+00
## Cumulative Proportion 0.3901 0.6847 0.8471 0.93195 1.00000 1.000e+00
ggplot(pca_df, aes(x = PC1, y = PC2, color = label)) +
geom_point() +
geom_text(aes(label = name), position = position_nudge(y = -2), show.legend = F, size = 3) +
scale_color_manual(values = colors_default) +
scale_x_continuous(expand = c(0.2, 0))

Correlation heatmap
pheatmap(
cor(rld$matrix),
annotation_col = as.data.frame(colData(dds)) %>% select(label),
color = brewer.pal(8, 'YlOrRd')
)

Wald test results
# DE analysis using Wald test
dds_full <- DESeq(dds)
colData(dds_full)
## DataFrame with 6 rows and 26 columns
## Assay Type AvgSpotLen Bases BioProject BioSample Bytes Center Name Consent DATASTORE filetype DATASTORE provider DATASTORE region Experiment treatment GEO_Accession (exp) Instrument LibraryLayout LibrarySelection LibrarySource Organism Platform label ReleaseDate Sample Name source_name SRA Study sizeFactor
## <character> <numeric> <numeric> <character> <character> <numeric> <character> <character> <character> <character> <character> <character> <factor> <character> <character> <character> <character> <character> <character> <character> <factor> <POSIXct> <character> <character> <character> <numeric>
## SRR13535276 RNA-Seq 300 8225466000 PRJNA694971 SAMN17588686 3252113587 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943360 A GSM5043430 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA in space without gravity 2021-09-09 GSM5043430 C2C12 proliferating .. SRP303354 0.970043
## SRR13535278 RNA-Seq 300 9203426700 PRJNA694971 SAMN17588684 3619152333 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943362 A GSM5043433 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA in space without gravity 2021-09-09 GSM5043433 C2C12 proliferating .. SRP303354 1.339835
## SRR13535280 RNA-Seq 300 9323939700 PRJNA694971 SAMN17588682 3735905901 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943364 A GSM5043436 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA in space without gravity 2021-09-09 GSM5043436 C2C12 proliferating .. SRP303354 1.088304
## SRR13535300 RNA-Seq 300 12820015200 PRJNA694971 SAMN17587361 5047533646 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943384 E GSM5043471 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA on land 2021-09-09 GSM5043471 C2C12 proliferating .. SRP303354 1.435633
## SRR13535302 RNA-Seq 300 12499917600 PRJNA694971 SAMN17587359 4941074444 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943386 E GSM5043475 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA on land 2021-09-09 GSM5043475 C2C12 proliferating .. SRP303354 0.769365
## SRR13535304 RNA-Seq 300 7150086300 PRJNA694971 SAMN17587357 2845819297 GEO public fastq,sra gs,ncbi,s3 gs.US,ncbi.public,s3.. SRX9943388 E GSM5043478 Illumina HiSeq 2500 PAIRED cDNA TRANSCRIPTOMIC Mus musculus ILLUMINA on land 2021-09-09 GSM5043478 C2C12 proliferating .. SRP303354 0.596305
# Wald test results
res <- results(
dds_full,
contrast = c('treatment', condition, control),
alpha = pval_cutoff
)
res
## log2 fold change (MLE): treatment A vs E
## Wald test p-value: treatment A vs E
## DataFrame with 14648 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000025900 4.91415 -4.920389399 2.012472 -2.44494841 0.0144873 NA
## ENSMUSG00000098104 4.09534 0.853116865 1.107035 0.77063255 0.4409248 NA
## ENSMUSG00000033845 107.62217 -0.107664037 0.423914 -0.25397623 0.7995139 0.928640
## ENSMUSG00000102275 2.36352 -0.391339523 1.464850 -0.26715328 0.7893511 NA
## ENSMUSG00000025903 97.37418 -0.000670711 0.485955 -0.00138019 0.9988988 0.999586
## ... ... ... ... ... ... ...
## ENSMUSG00000061654 1.69275 1.4242122 2.576074 0.552861 0.580358 NA
## ENSMUSG00000079834 28.80697 0.9987955 0.808449 1.235446 0.216665 0.549180
## ENSMUSG00000095041 184.20628 0.0979396 0.573811 0.170683 0.864473 0.954516
## ENSMUSG00000063897 31.54450 -0.2099739 0.718745 -0.292140 0.770180 0.915941
## ENSMUSG00000095742 10.11070 0.1499996 0.948601 0.158127 0.874357 0.957516
mcols(res)
## DataFrame with 6 rows and 2 columns
## type description
## <character> <character>
## baseMean intermediate mean of normalized c..
## log2FoldChange results log2 fold change (ML..
## lfcSE results standard error: trea..
## stat results Wald statistic: trea..
## pvalue results Wald test p-value: t..
## padj results BH adjusted p-values
summary(res)
##
## out of 14648 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 312, 2.1%
## LFC < 0 (down) : 184, 1.3%
## outliers [1] : 179, 1.2%
## low counts [2] : 2840, 19%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
plotDispEsts(dds_full)

Summary details
# Upregulated genes (LFC > 0)
res_sig_df %>% filter(log2FoldChange > 0)
# Downregulated genes (LFC < 0)
res_sig_df %>% filter(log2FoldChange < 0)
# Outliers (pvalue and padj are NA)
res[which(is.na(res$pvalue)), ]
## log2 fold change (MLE): treatment A vs E
## Wald test p-value: treatment A vs E
## DataFrame with 179 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000067780 318.8356 -2.370696 1.42870 -1.659339 NA NA
## ENSMUSG00000025981 152.0780 -0.293486 1.15037 -0.255123 NA NA
## ENSMUSG00000038349 100.7994 -3.942602 1.21130 -3.254840 NA NA
## ENSMUSG00000026024 50.6697 -3.517479 1.25525 -2.802213 NA NA
## ENSMUSG00000085842 21.9171 5.164547 2.11409 2.442915 NA NA
## ... ... ... ... ... ... ...
## ENSMUSG00000005871 403.6343 -0.524237 0.979146 -0.535403 NA NA
## ENSMUSG00000044595 41.1231 1.833181 1.681042 1.090503 NA NA
## ENSMUSG00000024597 346.2137 -1.294730 0.973882 -1.329453 NA NA
## ENSMUSG00000118138 23.2910 5.906357 3.529023 1.673652 NA NA
## ENSMUSG00000033417 290.9183 -0.928829 0.966631 -0.960893 NA NA
# Low counts (only padj is NA)
res[which(is.na(res$padj) & !is.na(res$pvalue)), ]
## log2 fold change (MLE): treatment A vs E
## Wald test p-value: treatment A vs E
## DataFrame with 2840 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000025900 4.91415 -4.920389 2.01247 -2.444948 0.0144873 NA
## ENSMUSG00000098104 4.09534 0.853117 1.10703 0.770633 0.4409248 NA
## ENSMUSG00000102275 2.36352 -0.391340 1.46485 -0.267153 0.7893511 NA
## ENSMUSG00000102135 4.94574 -0.328110 1.08090 -0.303554 0.7614681 NA
## ENSMUSG00000098201 2.01347 -0.896309 1.72861 -0.518514 0.6040997 NA
## ... ... ... ... ... ... ...
## ENSMUSG00000064344 2.73923 0.163918 1.41592 0.1157680 0.907836 NA
## ENSMUSG00000064349 3.00782 -0.127341 1.24975 -0.1018934 0.918841 NA
## ENSMUSG00000064358 2.70135 0.141492 1.65446 0.0855216 0.931847 NA
## ENSMUSG00000064369 4.23154 1.425316 1.24382 1.1459174 0.251829 NA
## ENSMUSG00000061654 1.69275 1.424212 2.57607 0.5528615 0.580358 NA
Shrunken LFC results
plotMA(res)

# Shrunken LFC results
res_shrunken <- lfcShrink(
dds_full,
coef = str_c('treatment_', condition, '_vs_', control),
type = 'apeglm'
)
res_shrunken
## log2 fold change (MAP): treatment A vs E
## Wald test p-value: treatment A vs E
## DataFrame with 14648 rows and 5 columns
## baseMean log2FoldChange lfcSE pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000025900 4.91415 -0.241445095 0.566057 0.0144873 0.152195
## ENSMUSG00000098104 4.09534 0.140171540 0.461172 0.4409248 0.746712
## ENSMUSG00000033845 107.62217 -0.059541700 0.318136 0.7995139 0.929527
## ENSMUSG00000102275 2.36352 -0.036424512 0.454994 0.7893511 NA
## ENSMUSG00000025903 97.37418 0.000553662 0.339890 0.9988988 0.999554
## ... ... ... ... ... ...
## ENSMUSG00000061654 1.69275 0.0466583 0.472044 0.580358 NA
## ENSMUSG00000079834 28.80697 0.2859022 0.497015 0.216665 0.552069
## ENSMUSG00000095041 184.20628 0.0410310 0.367468 0.864473 0.954845
## ENSMUSG00000063897 31.54450 -0.0639566 0.400050 0.770180 0.916739
## ENSMUSG00000095742 10.11070 0.0313236 0.425676 0.874357 0.958102
plotMA(res_shrunken)

mcols(res_shrunken)
## DataFrame with 5 rows and 2 columns
## type description
## <character> <character>
## baseMean intermediate mean of normalized c..
## log2FoldChange results log2 fold change (MA..
## lfcSE results posterior SD: treatm..
## pvalue results Wald test p-value: t..
## padj results BH adjusted p-values
summary(res_shrunken, alpha = pval_cutoff)
##
## out of 14648 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 308, 2.1%
## LFC < 0 (down) : 175, 1.2%
## outliers [1] : 179, 1.2%
## low counts [2] : 2272, 16%
## (mean count < 4)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
Summary details
# Upregulated genes (LFC > 0)
res_shrunken_sig_df %>% filter(log2FoldChange > 0)
# Downregulated genes (LFC < 0)
res_shrunken_sig_df %>% filter(log2FoldChange < 0)
# Outliers (pvalue and padj are NA)
res_shrunken[which(is.na(res_shrunken$pvalue)), ]
## log2 fold change (MAP): treatment A vs E
## Wald test p-value: treatment A vs E
## DataFrame with 179 rows and 5 columns
## baseMean log2FoldChange lfcSE pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000067780 318.8356 -0.2241170 0.536162 NA NA
## ENSMUSG00000025981 152.0780 -0.0416596 0.442087 NA NA
## ENSMUSG00000038349 100.7994 -2.9602838 1.538360 NA NA
## ENSMUSG00000026024 50.6697 -0.5579450 1.140736 NA NA
## ENSMUSG00000085842 21.9171 0.1717354 0.521148 NA NA
## ... ... ... ... ... ...
## ENSMUSG00000005871 403.6343 -0.0999653 0.440350 NA NA
## ENSMUSG00000044595 41.1231 0.1293951 0.487081 NA NA
## ENSMUSG00000024597 346.2137 -0.2710043 0.521619 NA NA
## ENSMUSG00000118138 23.2910 0.0577704 0.480737 NA NA
## ENSMUSG00000033417 290.9183 -0.1876146 0.470764 NA NA
# Low counts (only padj is NA)
res_shrunken[which(is.na(res_shrunken$padj) & !is.na(res_shrunken$pvalue)), ]
## log2 fold change (MAP): treatment A vs E
## Wald test p-value: treatment A vs E
## DataFrame with 2272 rows and 5 columns
## baseMean log2FoldChange lfcSE pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000102275 2.36352 -0.0364245 0.454994 0.789351 NA
## ENSMUSG00000098201 2.01347 -0.0631973 0.464941 0.604100 NA
## ENSMUSG00000103903 3.35339 0.0953761 0.477444 0.365967 NA
## ENSMUSG00000079671 1.77247 -0.0315545 0.459385 0.799751 NA
## ENSMUSG00000083422 2.09977 -0.1183241 0.484283 0.290668 NA
## ... ... ... ... ... ...
## ENSMUSG00000064342 2.42848 0.0625178 0.469763 0.546405 NA
## ENSMUSG00000064344 2.73923 0.0168427 0.451690 0.907836 NA
## ENSMUSG00000064349 3.00782 -0.0158357 0.444957 0.918841 NA
## ENSMUSG00000064358 2.70135 0.0107949 0.458129 0.931847 NA
## ENSMUSG00000061654 1.69275 0.0466583 0.472044 0.580358 NA
Visualizing results
Heatmaps
# Plot normalized counts (z-scores)
pheatmap(counts_sig_norm[2:7],
color = brewer.pal(8, 'YlOrRd'),
cluster_rows = T,
show_rownames = F,
annotation_col = as.data.frame(colData(dds)) %>% select(label),
border_color = NA,
fontsize = 10,
scale = 'row',
fontsize_row = 10,
height = 20)

# Plot log-transformed counts
pheatmap(counts_sig_log[2:7],
color = rev(brewer.pal(8, 'RdYlBu')),
cluster_rows = T,
show_rownames = F,
annotation_col = as.data.frame(colData(dds)) %>% select(label),
border_color = NA,
fontsize = 10,
fontsize_row = 10,
height = 20)

# Plot log-transformed counts (top 24 DE genes)
pheatmap(counts_sig_log %>% filter(ensembl_gene_id %in% (res_sig_df %>% head(24))$ensembl_gene_id) %>% select(-ensembl_gene_id) %>% column_to_rownames(var = 'mgi_symbol'),
color = rev(brewer.pal(8, 'RdYlBu')),
cluster_rows = T,
show_rownames = T,
annotation_col = as.data.frame(colData(dds)) %>% select(label),
fontsize = 10,
fontsize_row = 10,
height = 20)

Volcano plots
# Unshrunken LFC
res_df %>%
mutate(
sig_threshold = if_else(
padj < pval_cutoff & abs(log2FoldChange) >= lfc_cutoff,
if_else(log2FoldChange > 0, 'DE-up', 'DE-down'),
'non-DE'
)
) %>%
filter(!is.na(sig_threshold)) %>%
ggplot() +
geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = sig_threshold)) +
scale_color_manual(values = c('blue', 'red', 'gray')) +
xlab('log2 fold change') +
ylab('-log10 adjusted p-value')

# Shrunken LFC
res_shrunken_df %>%
mutate(
sig_threshold = if_else(
padj < pval_cutoff & abs(log2FoldChange) >= lfc_cutoff,
if_else(log2FoldChange > 0, 'DE-up', 'DE-down'),
'non-DE'
)
) %>%
filter(!is.na(sig_threshold)) %>%
ggplot() +
geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = sig_threshold)) +
scale_color_manual(values = c('blue', 'red', 'gray')) +
xlab('log2 fold change') +
ylab('-log10 adjusted p-value')

GSEA (all)
Hallmark genesets
# Shrunken LFC
get_fgsea_res(rank_lfc, mm_h) %>% plot_enrichment_table(rank_lfc, mm_h)

# Wald stat
get_fgsea_res(rank_stat, mm_h) %>% plot_enrichment_table(rank_stat, mm_h)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_h) %>% plot_enrichment_table(rank_pval, mm_h)

GO biological process
# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_bp) %>% plot_enrichment_table(rank_lfc, mm_c5_bp)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_bp) %>% plot_enrichment_table(rank_stat, mm_c5_bp)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_bp) %>% plot_enrichment_table(rank_pval, mm_c5_bp)

GO cellular component
# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_cc) %>% plot_enrichment_table(rank_lfc, mm_c5_cc)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_cc) %>% plot_enrichment_table(rank_stat, mm_c5_cc)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_cc) %>% plot_enrichment_table(rank_pval, mm_c5_cc)

GO molecular function
# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_mf) %>% plot_enrichment_table(rank_lfc, mm_c5_mf)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_mf) %>% plot_enrichment_table(rank_stat, mm_c5_mf)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_mf) %>% plot_enrichment_table(rank_pval, mm_c5_mf)

GSEA (DE)
Hallmark genesets
# Shrunken LFC
get_fgsea_res(rank_lfc, mm_h) %>% plot_enrichment_table(rank_lfc, mm_h)

# Wald stat
get_fgsea_res(rank_stat, mm_h) %>% plot_enrichment_table(rank_stat, mm_h)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_h) %>% plot_enrichment_table(rank_pval, mm_h)

GO biological process
# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_bp) %>% plot_enrichment_table(rank_lfc, mm_c5_bp)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_bp) %>% plot_enrichment_table(rank_stat, mm_c5_bp)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_bp) %>% plot_enrichment_table(rank_pval, mm_c5_bp)

GO cellular component
# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_cc) %>% plot_enrichment_table(rank_lfc, mm_c5_cc)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_cc) %>% plot_enrichment_table(rank_stat, mm_c5_cc)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_cc) %>% plot_enrichment_table(rank_pval, mm_c5_cc)

GO molecular function
# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_mf) %>% plot_enrichment_table(rank_lfc, mm_c5_mf)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_mf) %>% plot_enrichment_table(rank_stat, mm_c5_mf)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_mf) %>% plot_enrichment_table(rank_pval, mm_c5_mf)

System info
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /home/chan/mRNA_seq_pipeline/.snakemake/conda/9a19315a020c824d12f8055f7c009b0f/lib/libopenblasp-r0.3.18.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] fgsea_1.20.0 RColorBrewer_1.1-2 pheatmap_1.0.12 DESeq2_1.34.0 SummarizedExperiment_1.24.0 Biobase_2.54.0 MatrixGenerics_1.6.0 matrixStats_0.61.0 GenomicRanges_1.46.0 GenomeInfoDb_1.30.0 IRanges_2.28.0 S4Vectors_0.32.0 BiocGenerics_0.40.0 scales_1.1.1 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4 readr_2.1.1 tidyr_1.1.4 tibble_3.1.6 ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-2 ellipsis_0.3.2 XVector_0.34.0 fs_1.5.1 rstudioapi_0.13 farver_2.1.0 bit64_4.0.5 mvtnorm_1.1-3 AnnotationDbi_1.56.1 fansi_0.4.2 apeglm_1.16.0 lubridate_1.8.0 xml2_1.3.3 splines_4.1.0 cachem_1.0.6 geneplotter_1.72.0 knitr_1.35 jsonlite_1.7.2 broom_0.7.10 annotate_1.72.0 dbplyr_2.1.1 png_0.1-7 compiler_4.1.0 httr_1.4.2 backports_1.4.0 assertthat_0.2.1 Matrix_1.3-4 fastmap_1.1.0 cli_3.1.0 htmltools_0.5.2 tools_4.1.0 coda_0.19-4 gtable_0.3.0 glue_1.5.1 GenomeInfoDbData_1.2.7 fastmatch_1.1-3 Rcpp_1.0.7 bbmle_1.0.24 cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.3.8 Biostrings_2.62.0 xfun_0.28 rvest_1.0.2 lifecycle_1.0.1 XML_3.99-0.8 MASS_7.3-54 zlibbioc_1.40.0 vroom_1.5.7 hms_1.1.1 parallel_4.1.0 yaml_2.2.1 memoise_2.0.1 gridExtra_2.3 emdbook_1.3.12 bdsmatrix_1.3-4 stringi_1.7.6 RSQLite_2.2.8 highr_0.9 genefilter_1.76.0 BiocParallel_1.28.0 rlang_0.4.12 pkgconfig_2.0.3 bitops_1.0-7 evaluate_0.14 lattice_0.20-45 labeling_0.4.2 bit_4.0.4 tidyselect_1.1.1 plyr_1.8.6 magrittr_2.0.1 R6_2.5.1 generics_0.1.1 DelayedArray_0.20.0 DBI_1.1.1 pillar_1.6.4 haven_2.4.3 withr_2.4.3 survival_3.2-13 KEGGREST_1.34.0 RCurl_1.98-1.5 modelr_0.1.8 crayon_1.4.2 utf8_1.2.2 tzdb_0.2.0 rmarkdown_2.11 locfit_1.5-9.4 grid_4.1.0 readxl_1.3.1 data.table_1.14.2 blob_1.2.2 reprex_2.0.1 digest_0.6.29 xtable_1.8-4 numDeriv_2016.8-1.1 munsell_0.5.0